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with arbitrary base algorithms. 
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I. INTRODUCTION 

Because of an increasing demand in high-energy astrophysics, numerical modeling of rel¬ 
ativistic kinetic plasmas has been growing in importance. To date, many simulations on 
relativistic kinetic processes have been performed, such as the Rankine-Hugoniot problem 
across a relativistic shoclP, magnetic reconnection and kinetic instabilities in a relativisti- 
cally hot current sheet and the kinetic Kelvin-Helmholtz instability in a relativistic flow 
sheai^. In these simulations, one has to carefully set up ultrarelativistic bulk flows and/or 
relativistically hot plasmas in their rest frame. Loading velocity distribution function, i.e., 
initializing particle velocities by using random variables according to a relativistic distribu¬ 
tion function, is essentially important. 

In nonrelativistic particle simulations, it is quite natural to begin with a Maxwell- 
Boltzmann distribution (Maxwellian in short). To load the Maxwellian, the Box-Muller 
algorithm is widely used.*^ One can easily initialize a distribution with a bulk drift velocity, 
by applying an offset to the particle velocities. 

In relativistic simulations, it is natural to begin with a relativistic Maxwellian, also known 
as the Jiitter-Synge distribution function.l^*!^ In order to load it, perhaps the Sobol^ algo¬ 
rithm is the most popular, at least in Monte-Carlo simulation community. The algorithm 
was formally proposed by Sobol^^ in a Russian proceeding. Its key results are outlined in 
Pozdnyakov et al.^^^^. Meanwhile, it is not so clear how to initialize particles according to 
the relativistic shifted-Maxwellian or moving population of other distributions. To the best 
of our knowledge, the algorithms for the Jiitter-Synge distribution have not been applied 
to the relativistic shifted-Maxwellian. Several alternative algorithms have been proposed. 
Swisdak^ applied a rejection method for a log-concave distribution function. Melzani et 
al.l^ utilized a numerical cumulative distribution function and cylindrical transformation. 

In this research note, we describe numerical methods to load relativistic Maxwellians 
in particle simulations. We first describe two base algorithms to load stationary relativis¬ 
tic Maxwellian, the inverse transform method and the Sobol method.^ Next we apply the 
Lorentz transformation to obtain the relativistic shifted-Maxwellian. Simple rejection meth¬ 
ods are proposed to deal with the spatial part of the Lorentz transformation. We validate 
the algorithms by test problems, followed by discussions. 
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II. STATIONARY RELATIVISTIC MAXWELLIAN 


We consider relativistic Maxwell distributions (Jiittner-Synge distributiorP^ in the fol¬ 
lowing form, 


f{u)d^u = 


N 


4:7rm‘^cTK2{mc^ /T) 


exp 


'ymc 




( 1 ) 


where u = 'yv is the spatial components of the 4-velocity, v is the velocity, 7 = [1 - 
is the Lorentz factor, m is the rest mass, c is the light speed, T is the temperature, 
and K 2 {x) is the modihed Bessel function of the second kind. The normalization constant 
is set such that the number density is N = f f{u)d^u. Hereafter we set m = 1 and c = 1 
for simplicity. We use uppercases for fluid quantities and lowercases for particle properties 
throughout the paper. 

To generate u, we consider the spherical transformation Uy, u^) = {u sin 6 cos 99 , u sin 6 sin 73 , u cos 6). 
Then Equation yields 


f{u)du 


N 


a/1 -|- 

T 



( 2 ) 


In the special case of A = 1, one can read this equation as a probability function with 
respect to u. 


We generate u whose distribution follows Equation by either the inverse transform 


method (Sec. HA) or the Sobol method (Sec. HB). We will describe these methods in the 
next subsections. 


After we obtain m, we generate w on a spherical surface 1^1 = m in the momentum space. 
Using uniform random variables Yi(0 < Ai < 1) and ^ 2(0 < X 2 < 1), we set u in the 
following way. 


Ux = U ( 2 X 1 - 1) 

< Uy = cos{2tiX2) ■ (3) 

Uz = 2u^/xi{l^^^^sm.{2^lX2) 

Then we obtain a relativistic Maxwellian which follows Equation [TJ 
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A. Inverse transform method 

We consider the cumulative distribution function F{u) with a practical upper bound 

^max; 

Fiu) = ( y ( J f{u)du^ 

/ /■“ \ / /■“max . _1 

~ Uo ( Jq • ( 4 ) 

In the nonrelativistic limit of T 1, Umax = 5uth is sufficient, where Vth = is the ther¬ 
mal speed. In the relativistically hot case of T > 1, Equation [^behaves like oc exp{—u/T)u'^ 
for M 3> 1. This decays slower than the nonrelativistic limit of oc exp[—(n/nth)^]'^^, and so 
we increase the upper bound to Umax = 20T. We usually prepare a numerical table of F{u) 
with 2000 or more grid points. Using a uniform random variable X 3 , we compute 

u = F-\Xs) (5) 

by referring and interpolating the table. 


B. Sobol method 


Let US consider the gamma distribution. Its probability function P{x) is given by 

1 


P{x; a, b) = 




(x > 0 ), 


( 6 ) 


6 “ Gamma(a 

where a and b are free parameters and Gamma(a;) is the Gamma function. The gamma 
distribution with an integer parameter a can be generated by multiple random variables 
Aj’s (0 < Xj < 1 ) in the following way,^ 


X 

b 


5^hiW. 


(7) 


2=1 


Sobol^ noticed that the right hand side of Equationj^is similar to the third-order Gamma 
distributions, 

P(m; 3,T) = exp (^-(■g) 

For a certain T, we initialize u by using three random variables (X 4 ... Xq), 


u 


- = - In A 4 - InXg - InXg = - lnX 4 X 5 X 6 . 


(9) 
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Comparing the exponential parts in Equations and we obtain a relativistic Maxwellian 
by the rejection method. By using another random variable Xj, we accept the particle if 

exp - - -j > X 7 . ( 10 ) 

Then we obtain u which is distributed by Equation Using Equation this criteria can 
be modihed to 

VT+V^ < - TlnXy) = -TlnX4X5X6X7. (11) 

This leads to a simple form of the Sobol’s criterion,®^ 

rj^ -u^> 1 , ( 12 ) 


where r] = —TlnX 4 X 5 X 6 X 7 . Note that r] and u share the same variables X 4 ,X 5 , and Xq. 
Make sure to avoid zero in X 4 ... X 7 , because In 0 is undehned. Once Equation[^ (or Eq. 10) 
is met, we continue to the next step of the spherical scattering (Eq. |^. 

Comparing the normalization factors in Equation with N = 1 and Equation we 
obtain the overall efficiency of the rejection method as a function of 


1 


^2(1/T). 


(13) 


Figure shows the acceptance efficiency of the Sobol method, as a function of T. The 
efficiency quickly decreases for T —)■ 0, while it approaches to 1 for T —)■ cx). 


III. RELATIVISTIC SHIFTED-MAXWELLIAN 

A. Lorentz transformation 

Next we discuss general properties of the Lorentz transformation for particle distributions. 
We consider the transformation between two frames, S and S'. We assume that particles 
are stationary in the reference frame S, and then we switch to a moving frame S' at the 4- 
velocity (F, —F/3,0, 0). Without losing generality, we consider the transformation in the +x 
direction. In S', we observe the particle distribution, boosted by the 4-velocity (F, +F j3, 0, 0). 
We denote the observed properties in S' by the prime sign ('). 

As the total particle number is conserved, we recognize 

f{x, u) d^x d^u = f\x', u') d^x' d^u . (14) 
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FIG. 1. Acceptance efficiency of the Sobol algorithm as a function of the temperature 


The 


black squares show numerical results in Section IV 


Here, d^x = dxdydz is the spatial volume element. Using the time element dt in the same 
frame, we consider the 4-dimensional volume element of dt and d^x that is moving at the 
4-vector of u. The 4-dimensional position [t, x, y, z) follows the Lorentz transformation, and 
so the 4-volume element vector {dt, dx, dy, dz) also follows the Lorentz transformation. Since 
the Jacobian of the Lorentz transformation matrix A is 1, the 4-volume dtdz’x is conserved, 
i.e., dtd^x = dt' dz’x'. Since we deal with the it-moving volume, the time element dt is 
related to the canonical time element dr in the following way, dt = 'jdr. We similarly see 
dt' = 'y'dr. Therefore we obtain 


'-yd^x = 'y'd^x'. 


( 15 ) 


This also indicates the length contraction for the volume. The transformation is slightly 
different for d^u, because u is constrained by — 7 ^ = — 1 . Without losing 
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generality, one can consider the Lorentz transformation by (P, —P/3, 0, 0) in the +x direction: 


7' = P7(l + / 3 n^), 
du'^ = T{dux + /3d7), 

du'y = duy, du^ = duz- 

Under the condition of 7^7 = u^dux, we obtain 

d^u d?u' 


(16) 

(17) 

(18) 


7 


7 


(19) 


From Equations. 15 and 19, we obtain d^xd^u = d^x' d^u'. This ensures 


f{x,u) = f{x',u'). 

We obtain a relativistic shifted-Maxwellian by simply translating Equation 


( 20 ) 


f{u) = f'{u') = 


N 


A7iTK2il/T) V T 


7 


N 


exp I — 


r(7' - /?<) 


( 21 ) 

( 22 ) 


4irrA'2(l/r) V T 

Since we know nice algorithms (Sec. II), we would like to initialize the particle momentum 
u in the S frame, and then translate it to the S' frame by the Lorentz transformation. 


u —)■ u'. This procedure contains the momentum-space transformation (Eq. 19). However, 


it does not take care of the spatial part of the transformation, d^x —)■ d^x' (Eq. 15). Using 
the F-frame quantities, the distribution in S' appears to the observer in the following way. 


f'{u')d^u' = /(u)^— 


(23) 


We recognize a volume transform factor (yVb); because the element volume in S is not 
identical to the element volume in S'. This issue is also addressed by Melzani et al.^^. 


^-=V{1 + I3vx). 
7 


(24) 


One can also interpret that the number density is reciprocal to the volume size cx {d^x)~^ 
(See also Eq. 15). Since both spacial and momentum transformation (Equations [l^ and 19) 
depends on u, the factor differs from particle to particle. This may sound tricky, but the 
above formula describes what the observer looks at. We obtain very different results without 
the volume transformation, as will be shown in Section IV. 
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In this line, we briefly outline relativistic fluid properties. We assume isotropic Maxwellian 
distribution. From Equation the number flux 4-vector yields 




(25) 


We see = {N',N'V') = N{T,Tf3). Equations 25 ensures that follows the Lorentz 


transformation, i.e., N'^ = A^N°‘, where A is the Lorentz tensor. 
Similarly, the stress-energy tensor yields, 

r'"' = j f{u)uv‘‘A, 


(26) 


Clearly it follows the Lorentz transformation, i.e., T'^'' = A^A^T"^. In this case, 

= f2(^ + P)-P, (27) 

= T‘^{8 + P)I3\ (28) 

where £ = j f{u)'yd^u = N{[K3{1/T)/K 2 {l/T)] — T} is the internal energy density and 
P = / f {u)ux{ux/'y)(Pu = NT is the pressure in the rest frame. 

B. Volume transform methods 


Here, we describe simple methods to deal with the volume transform factor (Eq. 24) 


It is impossible to deal with this by adjusting the cell size in PIC simulation, because the 
transformation differs from particle to particle. One can also change the weight of particles. 
However, we prefer not to do so, because the ratio of the heaviest particle to the lightest 
particle could be very large. 

We propose to adjust the particle number by a rejection method. Using a random variable 
Ag (0 < Ag < 1), we accept the particle if the following condition is met, 

1 1 




(29) 


The left hand side ranges from 0 to 1. If the condition is not met, then we re-initialize 
the particle momentum. The factor (1/2P) can be absorbed in the normalization constant, 
because we usually know the value of 2PA before loading particles. The expected value E[x] 
of the acceptance efficiency is 50%, 


E 


2P Vy 


= p + = 0 - 5 - 


(30) 
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Fluid rest frame (S) 


Observer frame (S’) 



Acceptance factor 



Momentum part of the transformation 

7 ' ^ r (7 + (5ux) 
u'x ^ + / 37 ) 

Spatial part of the transformation 


FIG. 2. Lorentz transformation of a relativistic hot plasma distribution. The bottom panel illus¬ 
trates the flipping method, which is responsible for the spatial part of the Lorentz transformation. 


If S is not the fluid rest frame, E[nj;] 7 ^ 0 and so the efficiency may vary. 

We can further improve the efficiency in a special case of a symmetric distribution. When 
f{ux) = f{—Ux), we multiply the acceptance factor by 2 , 

i(^)=(l + /3^t.J. (31) 

The (1/F) factor is absorbed in the total particle number. We take advantage of the fact that 
the second term of the right hand side is odd function of (or v^). When is negative, 
the acceptance factor ranges between 0 < (1 — \l3vx\) < 1- We reject the particles at the 
probability of \livx\- On the other hand, when fivx is positive, the factor ranges between 
1 < (1 + I^Vx) < 2. We accept all particles. In addition, we interpret that we need to add 
another set of particles at the probability of \l3vx\- If f{ux) = f{—Ux), the number of the 
rejected particles and the number of the particles to be added are equal. We simply reverse 
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the sign of Ux of the rejected particles, and then we add them to the positive-/3na; side. This 
logic is schematically illustrated in the bottom of Figure 

We summarize the algorithm in the following way. If the following condition is met for a 
random variable Xg, 

- f^Vx > Xg, (32) 

then we change Ux —t —Ux, before computing Here we combine the two conditions of 


—I3vx < 0 and —I3vx > Xg to one condition (Eq. 32). The acceptance efficiency is 100%. We 


call it the flipping method (Eq. 32) to distinguish it from the rejection method (Eq. 29). 


IV. TEST PROBLEMS 


In order to validate the algorithms, we carry out several test problems. We initialize 10® 
particles in all cases. The black squares in Figure [T] show the acceptance efficiency of the 
Sobol method, as a function of T. They are in excellent agreement with the expected curve 

(Eq.Q. 

We then compute relativistic shifted-Maxwellian by using the Sobol method and the flip¬ 


ping method (Eq. 32). We set T = 1 and we boost the particles by the bulk Lorentz factor 


F = (1,1.1,10) in the +Ux direction. Figure compares numerical results and analytic 
distributions in the moving frame S", integrated over u'y and u'^. All distributions are nor¬ 
malized by J f'Su' = TN. The following analytic solution is obtained by using a cylindrical 
transformation («(,, u'y, u'^) = («(,, u'^ cos 0, u'^ sin 0) in Equation 


22 


- 27 r 


/«) = 


f' {u')u'j_d(j)du'j_ 


iV(F^rT^ + T) 


exp 


F(v^rT^-0O 


(33) 


2F2A:2(1/T) V T 

The numerical results are in excellent agreement with the analytic solutions. The stationary 
Maxwellian looks OK. As F increases, the distribution is stretched in the +u'^ direction. 

we see /'«) oc < exp[-(F(^l -F u'^ - 0O/T)] ^ < exp[-«/2Fr)] 


From Equation 


33 


for —)■ oo. Therefore, the slope on the boosted side becomes extremely flat. For F = 1.1, 
the numerical results on the right side (uj, ~ 20) look slightly noisier than on the other 
side (m(, ~ —8). This is probably a unfair comparison, because the right slope has more 
gridpoints than the left slope in the low-density range. For F = 10, the distribution is highly 
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stretched in +u'^. Outside the hgure, it still remains f'{u'^)/VN 4 x 10“^ at u'^ = 100. 

It will be challenging to initialize such a distribution by a direct rejection method in the S' 
frame, because we have to extend the parameter domain 2P times longer in +u'^. This gives 
us another motivation to initialize particles in S and then boost it to the S' frame. 

Next, we compute several fluid quantities in the moving frame S'. After initializing 
the particles, we compute the flow vector and the stress-energy tensor Then 

we evaluate the average velocity N"" /N'^ = and the average energy flux = 

P/3(T -|- P)/N. The former is a direct indicator of the bulk motion, and the latter, the 
energy flux, plays a decisive role to the system evolution. The results are presented in 
Table We change two key parameters, the bulk Lorentz factor P = (1.1,10,10^) and the 
relativistic temperature T = (0.1,1,10). In the T = 0.1 case, we use the inverse transform 
method (Sec. II A), because the efficiency of the Sobol method falls to ~ 0.001. We also 
test the T = 10 case without the volume transformation. This incorrect case is denoted 
by the asterisk sign (*). In Table the first rows show the computed results. The second 
rows show the relative error to analytic solutions. As can be seen, the results appear to be 
accurate, except for the rightmost columns. 

Without the volume transformation, we see that the average bulk speed is inaccurate in 
Table 1^ This is crucial to initialize a relativistic current sheelP^, in which relativistically hot 
populations carry the electric current. The energy flux is significantly distorted, too. The 
average energy flux without the volume transformation is 

Since {£ + P)/£ —)■ 4/3 for T S> 1, we lose 25% of the energy flux, regardless of the 
bulk speed (3. We can similarly evaluate the average energy density without the volume 
transformation. It deviates from the right value by a factor of [1 -|- (f and therefore 
the error approaches 25% for P S> 1 and T 3> 1. 

V. DISCUSSION AND SUMMARY 

We first reviewed two algorithms to initialize the stationary relativistic Maxwellian. In 
addition to the simple inverse transform method, we have formally reviewed the Sobol 
algorithm. In our experience, the inverse transform method is faster than the Sobol method. 
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u'x 


FIG. 3. Distribution functions of Lorentz-boosted Maxwellians as a function of u'^. Nu¬ 

merical results are overplotted on the analytic curves (Eq. [33]). We set T = 1 for all cases. 


because it only requires 3 random variables. We don’t see any problems, as long as we 
prepare 10^-10^ grids in the table. The algorithm can deal with any spherically-symmetric 
distributions. On the other hand, the Sobol method has a strong mathematical background. 
It is very simple, and so we can easily avoid a bug. The method is certainly slower than 
the inverse transform method, because it uses 6 random variables. However, this will not 
be a big deal, because we use these algorithms for initialization. The only problem is that 
the Sobol method becomes extremely inefficient for the nonrelativistic limit of T 1. In 
such a case, we simply switch from the Sobol method to the inverse transform method or 
the Box-Muller method. Another promising option is the log-concave rejection method, 
described in Section H and Appendix in Swisdak^. The algorithm uses 4 random variables, 
its acceptance efficiency is ~90%, and it is nearly insensitive to T. 


After initializing the stationary Maxwellian, we apply the volume transformation before 
boosting the particle momentum. We have proposed the two algorithms, the rejection 
method (Eq. [2^ and the flipping method (Eq. 32). They require one more random variable. 
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TABLE I. Computed fluid quantities and relative errors. 


N'^ 

TN 

II 

o 

1— ‘ 

T=1.0 

T=10 

* 

o 

T—1 

II 

r=i.i 

0.416532 

0.416708 

0.416566 

0.288896 


1.6 X lO""^ 

2.6 X lO""^ 

7.7 X 10-5 

0.307 

r=io 

0.994989 

0.994996 

0.994958 

0.975918 


1.6 X 10“® 

8.9 X 10“® 

2.9 X 10-5 

0.0192 

r=ioo 

0.999950 

0.999950 

0.999950 

0.999658 


9.3 X 10"^° 

4.0 X 10"® 

1.3 X 10-^ 

2.9 X lO-"^ 

'Y'lOx 

TW 

T=0.1 

T=1.0 

T=10 

T=10* 

r=i.i 

0.580438 

2.00167 

18.3798 

13.7357 


2.9 X 10"^ 

5.5 X lO""^ 

1.4 X 10-5 

0.252 

r=io 

12.6063 

43.4805 

398.147 

298.838 


2.6 X 10"® 

1.1 X lO""^ 

8.5 X 10-“^ 

0.250 

r=ioo 

126.674 

437.062 

4001.75 

3003.43 


1.5 X 10“^ 

9.2 X 10-5 

7.4 X 10-^ 

0.250 


The flipping method is our hrst choice. Since it accepts all particles, the overall efficiency 
is the same as the base algorithm for the stationary one. As a representative case, the 
Sobol method with the flipping method are summarized in the pseudocode in Table We 
emphasize that our volume transform methods are generic. The flipping method can be 
combined with power-law, waterbag, or any other distributions, as long as it is symmetric 
in Ux in the S frame. Even when the distribution is non-symmetric, we can divert to the 


rejection method (Eq. 29). The acceptance efficiency is typically 50%, but it works in any 
cases. Swisdak^ used the log-concave rejection method twice for the shifted Maxwellian. 
According to his article, the overall efficiency is ~ 80% insensitive to T. This is a very good 
result. However, his algorithm is specialized for Maxwellians or possibly other exponential- 
type distributions. In contrast, our simple methods can deal with any kind of distributions. 


Using the test problems, we have demonstrated that the combinations of the base methods 
and the flipping method excellently work. Without the volume transformation, we recognize 
signihcant errors up to 25% in the average energy flux. This is because the volume transform 
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TABLE II. Sobol algorithm with the flipping method. 


repeat 

generate Xi, X 2 , X 3 , X 4 , uniform on (0, 1] 

-TinX 1 X 2 X 3 

T] i -rinAiA2A3A4 

until Iff — v? > 

generate X 5 , Xq, Xy, uniform on [0, 1] 

Ux U (2A5 1) 

Uy 2 uy^A 5(1 — A 5 ) cos(27rA6) 

Uz <— 2uY^A5^1AAY^sin(27rA6) 

if {-^Vx > Xj), Ux i - Ux 

Ux ■<— r(Ua; + /3\/l + "U^) 
return Ux,Uy, Uz 


factor (Eq. 24) is no longer constant for P > 1 and T !:§> 1. 

In summary, we have described numerical algorithms to load relativistic Maxwellians in 
particle simulations. The inverse transform method and the Sobol method are useful to load 


the stationary Maxwellian. For shifted Maxwellian, the rejection method (Eq. 29) and the 
flipping method (Eq. [3^ take care of the spatial part of the Lorentz transformation. These 


methods are simple and physically-transparent. They can be combined with arbitrary base 
algorithms. We hope that these algorithms are useful in relativistic kinetic simulations in 
high-energy astrophysics. 
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